Summary

We simulated a trait controlled by 1,000 QTL 100 times. For each, one population was selected for an increase in the trait and one population was not selected at all (drift only). Selection or drift persisted for a variable numebr of generations, with a variable number of post-selection generations simulated. Herein, we use these simulations to test the power of our test statistic to identify selection on the trait.

Set preliminary values and load function.

These values will store results of the loop

Pvals.20.5 <- c()
source("Ghat.R")

Parameters: 20-selected, 5-not-selected

Simulations were performed on the IRCF Biocluster using the software QMSim. Each iteration of the loop will load data, manipulate data, conduct a test for selection, and record the p-value.

for(sim in 1:100){
    iter <- sprintf("%03i",sim)
    directory="/home/beissingert/data/PolySelect_bigsims/Intermittent/r_intermittent"

    #load genotypes and phenotypes
        pheno <- read.table(paste(directory, "/20Select_5Not_data_",iter,".txt",sep=""),header=T,stringsAsFactors=F)
    map <- read.table(paste(directory, "/lm_mrk_",iter,".txt",sep=""),header=T,stringsAsFactors=F)
    geno <- read.table(paste(directory, "/20Select_5Not_mrk_",iter,".txt-head",sep=""),header=F,stringsAsFactors=F,skip=1,sep="", colClasses=c("numeric","character")) # read genos as characters #
    qtlMap <- read.table(paste(directory, "/lm_qtl_",iter,".txt",sep=""),header=T,stringsAsFactors=F)

    #load allele frequencies
        freqs0<-read.table(paste(directory, "/20GenSelect_freq_mrk_",iter,".txt",sep=""),header=T,nrows=nrow(map),fill=T,stringsAsFactors=F) #NOTE POP TAKEN FROM
    freqs20<-read.table(paste(directory, "/20Select_5Not_freq_mrk_",iter,".txt",sep=""),header=F,skip={nrow(map)+1},fill=T,stringsAsFactors=F) #

    ### Manipulate genotypes by coding to -1,0,1 and so that markers are columns, individuals are rows.
    gen <- matrix(NA,nrow=nrow(map),ncol=nrow(geno))
    gen <- as.data.frame(gen)
    names(gen) <- geno[,1]
    for(i in 1:1000){
      #print(i)
      tmp <- as.numeric(unlist(strsplit(geno[i,2],split="")))
      tmp[which(tmp == 0)] <- -1
      tmp[which(tmp == 3 | tmp ==4)] <- 0
      tmp[which(tmp==2)] <- 1
      gen[,i] <- tmp
     }
    gen<-t(gen)
    gc()

    ### Calculate allele frequencies
    #Gen 0
    names(freqs0)[4]<- "Allele1"
    names(freqs0)[5]<- "Allele2"
    freqs0$Allele2[which(substr(freqs0$Allele1,1,1)==2)] <- "2:1.00000" # put this in the spot for the second allele
    freqs0$Allele1[which(substr(freqs0$Allele1,1,1)==2)] <- "1:0.00000" 
    freqs0$Allele2[which(substr(freqs0$Allele1,3,3)==1)] <- "2:0.00000" ##
    freqs0$Allele1 <- as.numeric(substr(freqs0$Allele1,3,1000))
    freqs0$Allele2 <- as.numeric(substr(freqs0$Allele2,3,1000))

    #Gen 20
    names(freqs20)[4]<- "Allele1"
    names(freqs20)[5]<- "Allele2"
    freqs20$Allele2[which(substr(freqs20$Allele1,1,1)==2)] <- "2:1.00000" # put this in the spot for the second allele
    freqs20$Allele1[which(substr(freqs20$Allele1,1,1)==2)] <- "1:0.00000" 
    freqs20$Allele2[which(substr(freqs20$Allele1,3,3)==1)] <- "2:0.00000"  ##
    freqs20$Allele1 <- as.numeric(substr(freqs20$Allele1,3,1000))
    freqs20$Allele2 <- as.numeric(substr(freqs20$Allele2,3,1000))

    #Calculate change
    change2<-freqs20$Allele2-freqs0$Allele2

### Setup for function
        geno <- gen
        phen <- as.matrix(pheno[1:1000,10])
        rownames(phen)<-pheno[1:1000,1]
    change=change2
    perms <- 10000
    blockSize=100000/100

### Determine num_eff
        potentials <- c(40000,13333,6667,4000,2857,2222,1818,1538,1333,1176,1053,667,400,286,222)
        ldTab <- read.table(paste(directory,"/20Select_5Not_ld_decay_",iter,".txt",sep=""),skip=43,nrows=15)
    ld <- as.numeric(substr(ldTab[,2],1,6))
    Meff <- potentials[which(ld <= 0.03)[1]]
    if(is.na(Meff)) Meff <- potentials[15]  #
    print(Meff)
    

    # Run function
    test<- Ghat_func(geno=geno,phen=phen,change=change2,method = "scale", num_eff = Meff,  perms=1000,plot="Both")
    Pvals.20.5[as.numeric(iter)] <- test$p.val
}
## [1] 400
## Loading required package: rrBLUP
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 667
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 667
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 667
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 286
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

## [1] 400
## [1] "Doing Permutation Test"

Report the results

Below are some summary statistics about the results. Note also that I save the workspace.

hist(Pvals.20.5)

print(Pvals.20.5)
##   [1] 1.382589e-05 2.337717e-08 1.131123e-05 5.328571e-09 6.049600e-08
##   [6] 4.500957e-08 1.477593e-04 1.534103e-08 8.376485e-06 2.479922e-06
##  [11] 6.809595e-07 4.008208e-08 7.051275e-06 4.410218e-09 1.323629e-08
##  [16] 1.623001e-05 1.186751e-09 2.247514e-07 1.995555e-05 3.057348e-07
##  [21] 1.385703e-05 7.575213e-09 5.230086e-06 1.994869e-10 7.539975e-06
##  [26] 6.406634e-08 4.089458e-09 2.423449e-09 7.493862e-08 3.641873e-12
##  [31] 3.225783e-05 1.894860e-05 2.882265e-07 1.828686e-08 7.477186e-05
##  [36] 1.606267e-05 1.640856e-07 1.829604e-06 5.793842e-07 2.244819e-05
##  [41] 7.061032e-06 9.341524e-06 2.058416e-05 2.861684e-04 6.405364e-05
##  [46] 6.664228e-07 2.249467e-07 2.094596e-06 2.370846e-09 5.173802e-09
##  [51] 2.278162e-06 2.602849e-06 3.985291e-10 9.442805e-07 1.230474e-05
##  [56] 1.068580e-08 4.882379e-07 1.600186e-04 4.254955e-07 2.918771e-09
##  [61] 6.553425e-06 1.262617e-06 1.368520e-07 2.047641e-06 4.820260e-06
##  [66] 2.257799e-05 1.556989e-04 4.156825e-12 3.087850e-04 1.422297e-03
##  [71] 5.909762e-05 7.818795e-06 1.956025e-07 4.274564e-07 2.178949e-06
##  [76] 1.037955e-07 2.435861e-05 7.090692e-05 5.432631e-06 5.187017e-05
##  [81] 3.667993e-04 4.589323e-05 7.803839e-07 3.726315e-08 1.658026e-06
##  [86] 1.833542e-06 2.319005e-06 8.725873e-04 2.980138e-05 4.504325e-07
##  [91] 4.490898e-06 8.184397e-05 2.504850e-03 6.434933e-07 3.618112e-09
##  [96] 2.089820e-08 8.721400e-09 4.695622e-05 1.039105e-07 2.740904e-03
length(which(Pvals.20.5<=0.05))/length(Pvals.20.5)
## [1] 1
save.image("on20off5.RData")